﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace GeodeticCoordinates
{
    internal class Calculate
    {
        private int _n; // 投影方式 3度带还是6度带
        public int n { set { _n = value; } get { return _n; } }

        private int _n0; // 带号
        public int n0 { set { _n0 = value; } get { return _n0; } }
        public struct Xishu
        {
            private double _a0;
            public double a0 { set { _a0 = value; } get { return _a0; } }

            private double _a1;
            public double a1 { set { _a1 = value; } get { return _a1; } }

            private double _a2;
            public double a2 { set { _a2 = value; } get { return _a2; } }

            private double _a3;
            public double a3 { set { _a3 = value; } get { return _a3; } }

            private double _a4;
            public double a4 { set { _a4 = value; } get { return _a4; } }

            private double _a5;
            public double a5 { set { _a5 = value; } get { return _a5; } }

            private double _a6;
            public double a6 { set { _a6 = value; } get { return _a6; } }

            private double _a7;
            public double a7 { set { _a7 = value; } get { return _a7; } }
        }

        /// <summary>
        /// 将文本中的度转换为度分秒形式
        /// </summary>
        /// <param name="d">角度的度形式</param>
        /// <returns>度分秒形式</returns>
        public static string TranDegreeToDMs(double d)
        {
            int Degree = Convert.ToInt16(Math.Truncate(d));//度
            d = d - Degree;
            int M = Convert.ToInt16(Math.Truncate((d) * 60));//分
            int S = Convert.ToInt16(Math.Round((d * 60 - M) * 60));
            if (S == 60)
            {
                M = M + 1;
                S = 0;
            }
            if (M == 60)
            {
                M = 0;
                Degree = Degree + 1;
            }
            string rstr = Degree.ToString() + ".";
            if (M < 10)
            {
                rstr = rstr + "0" + M.ToString();
            }
            else
            {
                rstr = rstr + M.ToString();
            }
            if (S < 10)
            {
                rstr = rstr + "0" + S.ToString();
            }
            else
            {
                rstr = rstr + S.ToString();
            }
            return rstr;
        }

        /// <summary>
        /// 度分秒转换为弧度制
        /// </summary>
        /// <param name="dms">度分秒</param>
        /// <returns>弧度制</returns>
        public static double Dms2H(double dms)
        {
            double Dms = dms + 0.000000001;
            double Du = Math.Truncate(Dms);
            double temp = (Dms - Du) * 100;
            double Fen = Math.Truncate(temp);
            double Miao = (temp - Fen) * 100 - 0.00001;

            double Du_1 = Du + Fen / 60.0 + Miao / 3600.0;
            return Du_1 * Math.PI / 180.0;
        }

        /// <summary>
        /// 弧度转度分秒
        /// </summary>
        /// <param name="Hu">弧度制</param>
        /// <returns>度分秒</returns>
        public static double H2Dms(double rad)
        {
            double degree = rad * 180 / Math.PI;
            double deg = Math.Truncate(degree);
            double min = Math.Truncate((degree - deg) * 60);
            double sec = Math.Round(((degree - deg) * 60 - min) * 60, 4);
            double dms = deg + min / 100 + sec / 10000;
            return dms;
        }

        /// <summary>
        /// 计算系数 高斯投影正算
        /// </summary>
        /// <param name="_e1">离心率</param>
        /// <returns>xishu的结构体形式</returns>
        public Xishu CalXishu(double _e1)
        {
            Xishu xishu = new Xishu();
            double e1 = Math.Sqrt(_e1);
            xishu.a0 = 1 + Math.Pow(e1, 2) * 3 / 4 + Math.Pow(e1, 4) * 45 / 64 + Math.Pow(e1, 6) * 175 / 256 + Math.Pow(e1, 8) * 11025 / 16384 + Math.Pow(e1, 10) * 43659 / 65536
                    + Math.Pow(e1, 12) * 693693 / 1048576 + Math.Pow(e1, 14) * 2760615 / 41944304 + Math.Pow(e1, 16) * 703956835 / 1073741824;
            xishu.a1 = Math.Pow(e1, 2) * 3 / 4 + Math.Pow(e1, 4) * 15 / 16 + Math.Pow(e1, 6) * 525 / 512 + Math.Pow(e1, 8) * 2205 / 2048 + Math.Pow(e1, 10) * 72765 / 65536
                + Math.Pow(e1, 12) * 297297 / 262144 + Math.Pow(e1, 14) * 19324305 / 16777216;
            xishu.a2 = Math.Pow(e1, 4) * 15 / 64 + Math.Pow(e1, 6) * 105 / 256 + Math.Pow(e1, 8) * 2205 / 4096 + Math.Pow(e1, 10) * 10395 / 16384
                + Math.Pow(e1, 12) * 1486485 / 2097152 + Math.Pow(e1, 14) * 6441435 / 8388608;
            xishu.a3 = Math.Pow(e1, 6) * 35 / 512 + Math.Pow(e1, 8) * 315 / 2048 + Math.Pow(e1, 10) * 31185 / 131072
                + Math.Pow(e1, 12) * 165165 / 524288 + Math.Pow(e1, 14) * 6441435 / 16777216;
            xishu.a4 = Math.Pow(e1, 8) * 315 / 16348 + Math.Pow(e1, 10) * 3465 / 65536 + Math.Pow(e1, 12) * 99009 / 1048576 + Math.Pow(e1, 14) * 585585 / 4194304;
            xishu.a5 = Math.Pow(e1, 10) * 693 / 131072 + Math.Pow(e1, 12) * 9009 / 524288 + Math.Pow(e1, 14) * 585585 / 16777216;
            xishu.a6 = Math.Pow(e1, 12) * 3003 / 2097152 + Math.Pow(e1, 14) * 45045 / 8388608;
            xishu.a7 = Math.Pow(e1, 14) * 6435 / 16777216;

            return xishu;

        }

        /// <summary>
        /// 高斯正算 计算中央子午线经度
        /// </summary>
        /// <param name="L">某地的经度</param>
        /// <param name="n">6度带还是3度带的投影方式</param>
        /// <returns>中央子午线经度以及所在带号</returns>
        public double[] GetL0(double L, int n)
        {
            double[] re = new double[2];
            double L0 = 0; // 中央子午线
            double dh = 0; // 带号
            if (n == 3)
            {
                dh = (int)((L + 1.5) / 3);
                L0 = dh * 3;//3°带中央子午线经度
            }
            else if (n == 6)
            {
                dh = (int)((L + 6) / 6);
                L0 = dh * 6 - 3;//6°带中央子午线经度
            }
            re[0] = dh;
            re[1] = L0;

            return re;
        }

        /// <summary>
        /// 高斯反算 中央子午线经度
        /// </summary>
        /// <param name="_n">带号</param>
        /// <param name="_n0">6度带还是3度带的投影方式</param>
        /// <returns>中央子午线经度</returns>
        public double GetL0(int _n, int _n0)
        {
            double re = 0;
            if (_n0 == 3)
            {
                re = _n * 3;
            }
            else if (_n0 == 6)
            {
                re = _n * 6 - 3;
            }
            return re;
        }

        /// <summary>
        /// 高斯正算
        /// </summary>
        /// <param name="angle"></param>
        /// <param name="l0"></param>
        /// <param name="a"></param>
        /// <param name="e1"></param>
        /// <param name="e2"></param>
        /// <returns></returns>
        public double[] CalLonLatToGaosi(double angle, double l0, double a, double e1, double e2)
        {
            {
                // 纬度B从度分秒换成弧度
                //double DMS = Convert.ToDouble(TranDegreeToDMs(angle));

                double DMS1 = angle;
                double radius = Dms2H(DMS1);
                double B = radius;

                // double L0 = ((int)((l0 - 1.5) / 3)) * 3; // 3°带中央子午线
                double[] temp_L0 = GetL0(l0, _n);
                double L0 = temp_L0[1]; // 中央子午线 度分秒的形式
                //double L0 = 111; // 中央子午线 度分秒的形式
                double l = (Dms2H(l0) - Dms2H(L0)) * 180 / Math.PI; //经度为度形式
                l = l * 3600; //换成秒的形式

                double N = a / Math.Sqrt(1 - e1 * Math.Pow(Math.Sin(B), 2));
                double t = Math.Tan(B);
                double yita = Math.Sqrt(e1 / (1 - e1)) * Math.Cos(B);


                Xishu xs = CalXishu(e1);
                e1 = Math.Sqrt(e1);

                double X = a * (1 - e1 * e1) * (
                    xs.a0 * B - xs.a1 * Math.Sin(2 * B) / 2 + xs.a2 * Math.Sin(4 * B) / 4 - xs.a3 * Math.Sin(6 * B) / 6 + xs.a4 * Math.Sin(8 * B) / 8
                    - xs.a5 * Math.Sin(10 * B) / 10 + xs.a6 * Math.Sin(12 * B) / 12 - xs.a7 * Math.Sin(14 * B) / 14);

                double x = X + N * Math.Sin(B) * Math.Cos(B) * l * l / 206265 / 206265 / 2
                + N * Math.Sin(B) * Math.Pow(Math.Cos(B), 3) * (5 - t * t + 9 * yita * yita + 4 * Math.Pow(yita, 4)) * Math.Pow(l, 4) / (24 * Math.Pow(206265, 4))
                + N * Math.Sin(B) * Math.Pow(Math.Cos(B), 5) * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(l, 6) / (720 * Math.Pow(206265, 6));


                double y = 500000 + N * Math.Cos(B) * l / 206265 + N * Math.Pow(Math.Cos(B), 3) * (1 - t * t + yita * yita) * Math.Pow(l, 3) / (6 * Math.Pow(206265, 3))
                + N * Math.Pow(Math.Cos(B), 5) * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * yita * yita - 58 * Math.Pow(t * yita, 2)) * Math.Pow(l, 5) / (120 * Math.Pow(206265, 5));

                double[] xy = new double[2];
                xy[0] = x;
                xy[1] = y;

                return xy;
            }
        }
        public double[] CalGaosiToLonLat(double _x, double _y, double _a, double _e1)
        {
            double e1 = Math.Sqrt(_e1);
            double a0 = 1 + Math.Pow(e1, 2) * 3 / 4 + Math.Pow(e1, 4) * 45 / 64 + Math.Pow(e1, 6) * 175 / 256 + Math.Pow(e1, 8) * 11025 / 16384 + Math.Pow(e1, 10) * 43659 / 65536
                    + Math.Pow(e1, 12) * 693693 / 1048576 + Math.Pow(e1, 14) * 27606215 / 41944304 + Math.Pow(e1, 16) * 703956835 / 1073741824;
            double C1 = Math.Pow(e1, 2) * 3 / 8 + Math.Pow(e1, 4) * 3 / 16 + Math.Pow(e1, 6) * 213 / 2048 + Math.Pow(e1, 8) * 255 / 4096;
            double C2 = Math.Pow(e1, 4) * 21 / 256 + Math.Pow(e1, 6) * 21 / 256 + Math.Pow(e1, 8) * 533 / 8192;
            double C3 = Math.Pow(e1, 6) * 151 / 6144 + Math.Pow(e1, 8) * 151 / 4096;
            double C4 = Math.Pow(e1, 8) * 1097 / 131072;
            double a = _a;

            double x = _x;
            double y = _y - 500000;

            double B0 = x / (a * (1 - e1 * e1) * a0);
            double Bf = B0 + C1 * Math.Sin(2 * B0) + C2 * Math.Sin(4 * B0) + C3 * Math.Sin(6 * B0) + C4 * Math.Sin(8 * B0);

            double tf = Math.Tan(Bf);
            double yitaf = Math.Sqrt(e1 * e1 / (1 - e1 * e1)) * Math.Cos(Bf);
            double W = Math.Sqrt(1 - Math.Pow(e1 * Math.Sin(Bf), 2));
            double Mf = a * (1 - e1 * e1) / Math.Pow(W, 3);
            double Nf = a / Math.Sqrt(1 - Math.Pow(e1 * Math.Sin(Bf), 2));

            double L0 = GetL0(_n0, _n); //中央子午线

            double B = Bf - tf * Math.Pow(y, 2) / (24 * Mf * Nf) + Math.Pow(y, 4) * tf * (5 + 3 * tf * tf + yitaf * yitaf - 9 * tf * tf * yitaf * yitaf) / (24 * Mf * Math.Pow(Nf, 3))
                + Math.Pow(y, 6) * tf * (61 + 90 * tf * tf + 45 * Math.Pow(tf, 4)) / (720 * Mf * Math.Pow(Nf, 5));

            B = Convert.ToDouble(H2Dms(B));

            double l = y / (Nf * Math.Cos(Bf)) - Math.Pow(y, 3) * (1 + 2 * tf * tf + yitaf * yitaf) / (6 * Math.Pow(Nf, 3) * Math.Cos(Bf))
                + Math.Pow(y, 5) * (5 + 28 * tf * tf + 6 * yitaf * yitaf + 24 * Math.Pow(tf, 4) + 8 * tf * tf * yitaf * yitaf) / (120 * Math.Pow(Nf, 5) * Math.Cos(Bf)); //弧度
            double L = Convert.ToDouble(TranDegreeToDMs(L0 + l * 180 / Math.PI));

            double[] xy = new double[2];
            xy[0] = B;
            xy[1] = L;

            return xy;
        }

        //public void GaosiZsuanFsuan(double _x, double _y, double _a, double _e1, double angle, double l0, double e2)
        //{
        //    switch (_n)
        //    {
        //        case 3:
        //            CalGaosiToLonLat(_x, _y, _a, _e1); //反算
        //            CalLonLatToGaosi(angle, l0, _a, _e1, e2); //正算
        //            break;
        //        case 6:
        //            CalGaosiToLonLat(_x, _y, _a, _e1); //反算
        //            CalLonLatToGaosi(angle, l0, _a, _e1, e2); //正算
        //            break;
        //    }
        //}
    }
}
